Basic Overview of the Problem

Why are we doing this pairwise?

Why not do full, joint pedigree reconstruction?

After all, doing so let’s you use all the available data

Because:

1. CKMR works best in sparse situations
2. Little advantage to “joint”, multi-individual inference (Colony, MasterBayes, Pedfactory) when true relationships are few and far between.
3. (Joint inference just takes a lot longer…)
4. Much harder to assess False Positive Rates in joint-inference,
   full-pedigree reconstruction methods than in pairwise inference.

Getting Statistical about “Are you my mother?”

Cast the question as a statistical test applied to pairs:

• \(H_0\) the pair is Unrelated
• \(H_A\) the pair is Parent-offspring (or full-sib, etc.)

Then, if it were the case that all pairs were either Unrelated or Parent-offspring

And, if we fixed our false positive rate (FPR) to be \(\alpha\)…

Then, amongst all possible statistical tests, the one with the highest power (lowest false negative rate) will be based on the likelihood ratio. \[ \frac{P(G_i, G_j~|~i~\mbox{and}~j~\mbox{are}~\mathrm{PO})} {P(G_1, G_2~|~i~\mbox{and}~j~\mbox{are}~\mathrm{Unrelated})} \] Where \(G_1\) and \(G_2\) denote the observed genotypes in the two individuals of the pair.

This holds for any two relationships: \[ \frac{P(G_i, G_j~|~i~\mbox{and}~j~\mbox{are}~R_1)} {P(G_1, G_2~|~i~\mbox{and}~j~\mbox{are}~R_2)} \] is a most powerful test (Neyman-Pearson Lemma).

Review

We previously showed we can use \(\boldsymbol{\kappa}(R)\) (the \(\kappa\) coefficients for relationship \(R\)) to calculate the joint probability \(P(G_1, G_2~|~i~\mbox{and}~j~\mbox{are}~R)\):

\[ \begin{aligned} P(G_i, G_j|\boldsymbol{\kappa}) &= \kappa_0 P(G_i)P(G_j) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\tiny\mathrm{no~gene~copies~IBD}} \\ &+ \kappa_1 P(G_i)P(G_j | G_i=\mathrm{ma}, G_j=\mathrm{kid}) ~~~~~~~~~{\tiny\mathrm{1~gene~copy~IBD}} \\ &+ \kappa_2 P(G_i)\mathcal{I}\{P(G_j)=P(G_i)\} ~~~~~~~~~~~~~~~~~~~~~~~~~~{\tiny\mathrm{2~gene~copies~IBD}} \\ \end{aligned} \] where \(\mathcal{I}\{\cdot\}\) is an indicator function returning 1 when what is inside it is true, and 0 otherwise.

What about at multiple loci?

With \(L\) different loci, writing the genotype at the \(\ell^\mathrm{th}\) locus as \(G^{(\ell)}\) we can use: \[ P(\boldsymbol{G}_1, \boldsymbol{G}_2~|~\boldsymbol{\kappa}(R)) = \prod_{\ell = 1}^L P(G_{1}^{(\ell)}, G_{2}^{(\ell)}~|~\boldsymbol{\kappa}(R)) \] This is actually the correct probability so long as:

1. The markers are not in linkage disequilibrium
2. The markers are not physically linked to one another

• It is pretty easy to verify that (1) is the case—merely test for LD.

• Assessing (2) is more difficult.

  • We might not have any idea where the markers are.
  • Even if we did, the calculation accounting for physical linkage is much more complicated.

So, assumption (2) often gets made for this calculation, even when it is known to be wrong. (We come back to this later.)

Back to likelihood ratios

Using that (approximate) joint probability of the genotypes we can compute the (approximate) likelihood ratio suggested by the NP-lemma: \[ \mathrm{For~PO~vs~Unrelated~pairs:}~~~~~~~~~~~~~~~~~~~\frac{P(\boldsymbol{G}_i, \boldsymbol{G}_j~|~\boldsymbol{\kappa}= (0, 1,0))} {P(\boldsymbol{G}_i, \boldsymbol{G}_j~|~\boldsymbol{\kappa}= (1, 0,0))} \]

Let’s give the log of that a name: \(\Lambda\):

\[ \Lambda_{PO/U} = \log\biggl(\frac{P(\boldsymbol{G}_1, \boldsymbol{G}_2~|~\boldsymbol{\kappa}= (0, 1,0))} {P(\boldsymbol{G}_1, \boldsymbol{G}_2~|~\boldsymbol{\kappa}= (1, 0,0))} \biggr) \]

It will be convenient later to recognize this as a sum:

\[ \begin{aligned} \Lambda_{PO/U} & = & \log\biggl( \prod_{\ell = 1}^L \frac{P(G_{1}^{(\ell)}, G_{2}^{(\ell)}~|~\boldsymbol{\kappa}= (0, 1, 0))} {P(G_{1}^{(\ell)}, G_{2}^{(\ell)}~|~\boldsymbol{\kappa}= (1, 0, 0))} \biggr) \\ & = & \sum_{\ell = 1}^L\log\biggl( \frac{P(G_{1}^{(\ell)}, G_{2}^{(\ell)}~|~\boldsymbol{\kappa}= (0, 1, 0))} {P(G_{1}^{(\ell)}, G_{2}^{(\ell)}~|~\boldsymbol{\kappa}= (1, 0, 0))} \biggr) \end{aligned} \]

Wait! What about Genotyping error?

• Good point.

• In our survey of genotyping technologies, we noted that they all can produce errors to some extent.

• Errors are most problematic in assessing likelihoods for the PO relationship.

• These probabilities can be computed in a way that accounts for genotyping error: \[ \Lambda_{PO/U} = \log\biggl(\frac{P(\boldsymbol{G}_1, \boldsymbol{G}_2~|~\boldsymbol{\mu}, \boldsymbol{\kappa}= (0, 1,0))} {P(\boldsymbol{G}_1, \boldsymbol{G}_2~|~\boldsymbol{\mu}, \boldsymbol{\kappa}= (1, 0,0))} \biggr) \]

• Error model can be simple (or complex)

• R package ‘CKMRsim’ (we will use it this afternoon) can accommodate all possible single-locus genotyping models.

How do we use these \(\Lambda\)’s?

• Truly PO pairs tend to have higher values of \(\Lambda_{PO/U}\) than do unrelated pairs.
• The distribution of \(\Lambda\) is not analytically tractable
• So, investigate it through Monte Carlo simulation.

So, Set a threshold \(\Lambda_c\). For example:

FALSE NEGATIVE RATE:

probability that a true kin pair (PO in this case) has \(\Lambda \leq \Lambda_c\).

FALSE POSITIVE RATE:

probability that a non kin pair (U in this case) has \(\Lambda > \Lambda_c\).

Hopes and Desires

• Obviously you want FNR as low as possible, but FPR must be low as well.
• How low must FPR be? That depends on the number of pairwise comparisons.
• 10,000 candidate parents and 10,000 candidate offspring = \(10^{8}\) pairwise comparisons.
• I think you should aim for FPR at least 100 times smaller than the reciprocal of the number of pairwise comparisons; preferably smaller.
• How do you estimate such small probabilities?

Calculating the FPR

“Vanilla” Monte Carlo Formulation, with MC sample size \(M\):

\[ \mathrm{FPR}_{\Lambda_c} = \sum_{\Lambda_{PO/U}}\mathcal{I}\{ \Lambda_{PO/U} > \Lambda_c\} P(\Lambda_{PO/U}~|~i,j~\mathrm{Unrelated}) \approx \frac{1}{M}\sum_{m=1}^M \mathcal{I}\{ \Lambda_{PO/U}^{(m)} > \Lambda_c\} \] with each \(\Lambda^{(m)}_{PO/U}\) simulated from its distribution assuming the pair is unrelated.

Importance sampling formulation rewrites the original expectation:

\[ \begin{aligned} \mathrm{FPR}_{\Lambda_c} &= \sum_{\Lambda_{PO/U}}\mathcal{I}\{ \Lambda_{PO/U} > \Lambda_c\} \frac{P(\Lambda_{PO/U}~|~i,j~\mathrm{Unrelated})} {P(\Lambda_{PO/U}~|~i,j~\mathrm{PO})} P(\Lambda_{PO/U}~|~i,j~\mathrm{PO}) \\ & \approx \frac{1}{M}\sum_{m=1}^M \biggl[\mathcal{I}\{ \Lambda_{PO/U}^{(m)} > \Lambda_c\} \frac{P(\Lambda^{(m)}_{PO/U}~|~i,j~\mathrm{Unrelated})} {P(\Lambda^{(m)}_{PO/U}~|~i,j~\mathrm{PO})}\biggr] \end{aligned} \] with each \(\Lambda^{(m)}_{PO/U}\) simulated from the importance sampling distribution, \(P(\Lambda^{(m)}_{PO/U}~|~i,j~\mathrm{PO})\).

• Useful for calculating FPR for comparison of any relationship to Unrelated,
  regardless of physical linkage.

Decreasing the FPR—Three main ways

1. Increase \(\Lambda_c\), thus increasing the FNR, too.

Easy to accommodate this in CKMR models if you can estimate the FNR.

Decreasing the FPR—Three main ways

2. Increasing the number of loci. Example HS vs Unrelated:

Decreasing the FPR—Three main ways

3. Increasing the heterozygosity of each locus
   • Using markers with more alleles:

Decreasing the FPR—Three main ways

3. Increasing the heterozygosity of each locus
   • Or, just by having allele frequencies that are closer to even.
   • Example: 100 SNPs at minor allele frequency of 0.05, 0.1, 0.2, 0.3, 0.5:

Physical Linkage reduces power

• But how does it do that?

• First, a digression into HS, A-N, and GP-GC, pairs.

• Note that even if you don’t have a map, you can still compute \(\Lambda\). It is a valuable statistic even if it is not the correct probability.

But how does physical linkage affect the distribution of \(\Lambda\)?

• Up till now, all simulations have assumed no linkage.
• This means that each locus, independently, is simulated to have 0, 1, or 2 gene copies IBD, according to \(\boldsymbol{\kappa}\).
• In reality, the IBD status of nearby sites in the genome will be correlated, and so their contributions to the sum: \[ \sum_{\ell = 1}^L\log\biggl( \frac{P(x_{1,\ell}, x_{2,\ell}~|~\boldsymbol{\kappa}_\mathrm{numer}}{P(x_{1,\ell}, x_{2,\ell}~|~\boldsymbol{\kappa}_\mathrm{denom}}\biggr) \] will also be positively correlated.
• So the variance of that sum is increased. Think \(\mathrm{Var}(X + Y)\). So there will be more spread in \(\Lambda\).

EXCEPT FOR SOME INTERESTING CASES

• The effect of linkage occurs due to variance in the total fraction of genome IBD at 0, 1, or 2 gene copies.
• Three special cases have no variance in that fraction:
  • U: \(\boldsymbol{\kappa}= (1, 0, 0)\)
  • PO: \(\boldsymbol{\kappa}= (0, 1, 0)\)
  • MZ: \(\boldsymbol{\kappa}= (0, 0, 1)\) (monozygotic twin, or “self”) And so, they are completely unaffected by physical linkage.
• So, \(\Lambda_{PO/U}\) is not affected by physical linkage.

However, distributions of \(\Lambda_{FS/U}\) or \(\Lambda_{HS/U}\) are affected by linkage

• But only their distributions when the truth is not U, PO, or MZ.
• Witness \(\Lambda_{FS/U}\) at 165 microhaps:

Have scanned 10000 pairs in Mendel output file.
Have scanned 10000 pairs in Mendel output file.

Let’s compute the mean and variance of those:

# A tibble: 8 × 5
# Groups:   simmed_from, numer_wts [8]
  simmed_from numer_wts denom_wts mean_lograt var_lograt
  <chr>       <chr>     <chr>           <dbl>      <dbl>
1 FS-linked   FS=1      U=1              44.8      192. 
2 FS-unlinked FS=1      U=1              44.7      117. 
3 HS-linked   HS=1      U=1              12.3       40.5
4 HS-unlinked HS=1      U=1              12.2       29.8
5 U-linked    FS=1      U=1             -36.9       64.0
6 U-linked    HS=1      U=1             -10.3       18.1
7 U-unlinked  FS=1      U=1             -37.0       67.0
8 U-unlinked  HS=1      U=1             -10.3       18.7

And that extra variance is causes the power to go down:

For unlinked:

# A tibble: 5 × 10
    FNR      FPR       se num_no…¹ Lambd…² pstar mc_me…³ numer…⁴ denom…⁵ true_…⁶
  <dbl>    <dbl>    <dbl>    <int>   <dbl> <chr> <chr>   <chr>   <chr>   <chr>  
1  0.01 3.69e-12 3.96e-13     9900    20.7 FS    IS      FS      U       U      
2  0.05 1.53e-14 9.52e-16     9500    27.3 FS    IS      FS      U       U      
3  0.1  5.74e-16 2.76e-17     9000    31.2 FS    IS      FS      U       U      
4  0.2  1.10e-17 4.26e-19     8000    35.6 FS    IS      FS      U       U      
5  0.3  4.87e-19 1.83e-20     7000    38.7 FS    IS      FS      U       U      
# … with abbreviated variable names ¹​num_nonzero_wts, ²​Lambda_star, ³​mc_method,
#   ⁴​numerator, ⁵​denominator, ⁶​true_relat

For linked:

# A tibble: 5 × 10
    FNR      FPR       se num_no…¹ Lambd…² pstar mc_me…³ numer…⁴ denom…⁵ true_…⁶
  <dbl>    <dbl>    <dbl>    <int>   <dbl> <chr> <chr>   <chr>   <chr>   <chr>  
1  0.01 1.16e- 9 3.15e-10     9989    12.9 FS    IS      FS      U       U      
2  0.05 1.16e-12 1.08e-13     9850    22.3 FS    IS      FS      U       U      
3  0.1  1.92e-14 1.23e-15     9525    27.0 FS    IS      FS      U       U      
4  0.2  1.09e-16 4.69e-18     8639    33.0 FS    IS      FS      U       U      
5  0.3  2.01e-18 7.59e-20     7465    37.3 FS    IS      FS      U       U      
# … with abbreviated variable names ¹​num_nonzero_wts, ²​Lambda_star, ³​mc_method,
#   ⁴​numerator, ⁵​denominator, ⁶​true_relat